Gutzwiller scheme for electrons and phonons: the half-filled Hubbard-Holstein model 
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We analyze the ground-state properties of strongly-correlated electrons coupled with phonons by 
means of a generalized Gutzwiller wavefunction which includes phononic degrees of freedom. We 
study in detail the paramagnetic half-filled Hubbard-Holstein model, where the electron-electron in- 
teraction can lead to a Mott transition, and the electron-phonon coupling to a bipolaronic transition. 
We critically discuss the quality of the proposed wavefunction in describing the various transitions 
and crossovers that occur as a function of the parameters. Previous variational attempts are recov- 
ered as particular choices of the wavefunction, while keeping all the variational freedom allows to 
access regions of the phase diagram otherwise inaccessible within previous variational approaches. 



I. INTRODUCTION 

Understanding the interplay between electron-phonon 
(e-ph) and strong electron-electron (e-e) interactions can 
be important to understand the properties of many com- 
pounds, such as manganites^ or fullerides 2 -. This problem 
has recently received a renewed attention, after a number 
of experiments on high-T c superconducting cuprates tes- 
tified for a non trivial role of the e-ph coupling in the un- 
derdoped region, where correlation effects are important 
and standard e-ph theories are hardly reliable 3 -^^. In- 
deed a full theoretical understanding of lattice effects in 
strongly-correlated systems is still lacking when the two 
interactions have intermediate or large strength. 

The essential features of a strongly-correlated system 
are captured in the celebrated Hubbard model 8 - and in 
its extensions. On the other hand, polaron formation 
due to e-ph coupling has been extensively studied in the 
absence of e-e repulsion, mainly in the limit of low elec- 
tron density within the molecular-crystal model first in- 
troduced by Holstein 9 . Nonperturbative approaches such 
as Quantum Monte Carlo, Density Matrix Renormaliza- 
tion Group and Dynamical Mean-Field Theory (DMFT) 
have enormously contributed to the understanding of 
these models and have been successfully applied to 
the Hubbard-Holstein mode&ALl^J^d^.. Nonetheless 
these methods require a considerable numerical effort, 
and it is not unfair to say that approximate analytic re- 
sults, even if less accurate, are still highly desirable. 

From this point of view, variational approaches are 
a very natural choice as they provide immediate in- 
formation on the ground-state properties. In the field 
of strongly-correlated systems, the trial wavefunction 
introduced by Gutzwiller 1 - and then generalized by 
Biinemann and coworkers^ 7 - has proved quite accurate 
for electronic systems with local interactions. The un- 
derlying recipe is to construct the correlated wavefunc- 
tion starting from an adequate uncorrelated one and 
projecting out the energetically unfavored states. The 
equivalence with the Kotliar-Ruckenstein (KR) slave- 



boson saddle-point solution^ 8 - allows also to go beyond 
the mean-field approximation in a controlled way. For 
the Holstein model several trial wavefunctions have been 
proposed generalizing the Lang-Firsov transformation^ 9 -, 
which yields the exact solution in the atomic limit, 
and it is believed to capture the essence of polaronic 
physics, i.e., the entanglement between electronic and 
phonic degrees of freedom. These kind of approaches 
have been successfully applied to the Holstein model, 
mainly in the limit of low carrier density or for spin- 
less electron s 20 ' 21 ' 22 ' 23 ' 24 ' 25 . In fact, a further complica- 
tion arises at finite density, since the e-ph coupling in- 
duces an effective c-e attractive interaction that cannot 
be handled in an independent particle picture. When the 
Hubbard repulsion is included, it directly competes with 
such phonon-mediated attraction. 

The typical strategy to treat correlated fermions cou- 
pled with phonons consists in deriving an effective model 
for electrons by averaging the full Hamiltonian over 
a suitable trial phonon wavefunction (which also im- 
plies the neglect of any electron-multiphonon residual 
interaction)™, and then to resort to some appropriate 
technique in order to tackle the electronic effective in- 
teraction. This can be done for instance by means of a 
Hartree-Fock decoupling scheme 2 ^ or exploiting the more 
accurate KR slave-boson mean fiel d 27 ! 28 ' 29 suitable for 
local e-e interactions. As discussed in details in Refs. 
I24II30I the assumption of a factorized wavefunction be- 
comes questionable in the regime where retardation ef- 
fects between the motion of electrons and lattice defor- 
mations imply a strong exchange of momentum. 

In order to better describe the interplay between lattice 
vibrations and local electronic configurations, we intro- 
duced in a previous paper 31 a Gutzwiller Phonon Wave- 
function (GPW) to account for the e-ph coupling be- 
sides strong electronic correlation and benchmarked it 
in the infinite-?/ limit of the Hubbard-Holstein model 
by comparing with exact DMFT calculations. In this 
limit a finite phonon-driven attraction is obviously inef- 
fective against the infinite repulsion, and the only effect 
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of phonons is the formation of polarons. 

In this paper we turn our attention to the finite-C/ 
regime, analyzing in details the paramagnetic solution 
(i.e., neglecting broken-symmetry phases). Relaxing the 
infinite-?/ assumption, the competition between repul- 
sion and attraction becomes effective, leading to a rich 
phase diagram. We focus on the half-filling case, in which 
the e-e and e-ph interactions can drive, respectively, the 
Mott and bipolaronic metal-insulator transition, so that 
their interplay is both more transparent and more qual- 
itative than for arbitrary fillings. Out of half- filling the 
system is indeed always metallic in the absence of sym- 
metry breaking, and the effect of the competition is es- 
sentially quantitative and therefore less clear—. 

The key property of our GPW wavefunction is to treat 
simultaneously electrons and phonons from the onset 
by introducing a generalized Gutzwiller projector where 
phonon quantum states are determined by the local elec- 
tronic configurations rather than by the average electron 
density. This allows to address the modification induced 
by the electronic correlation onto the lattice groundstate. 
Therefore, even if the multiphonon adiabatic regime of 
the weakly-correlated system is probably still misrepre- 
sented, we expect to capture the essential features of 
the model in the presence of a sizeable repulsive inter- 
action, when the relevant physics is essentially local be- 
cause the electron motion is slowed down by correlations. 
We derive a set of variational equations for the phonon 
wavefunctions and show that standard approaches can 
be obtained as particular choices for the solution of such 
equations. We obtain a phase diagram which strongly de- 
pends on the adiabaticity regime, and critically discuss 
our results and the quality of the proposed wavefunction. 

The paper is organized as follows. In Section |TT] we 
introduce the model Hamiltonian and the relevant equa- 
tions for the general case of arbitrary electron filling, and 
then we specialize them for the half-filling case we are in- 
terested in. At the end of this section we derive standard 
variational approaches for the Holstein problem as partic- 
ular variational ansatzs on the Gutzwiller wavefunction. 
In Section HTTl we present the general solution provided by 
our method for the half-filled system, whose properties 
are analyzed in the following Section HVl We then discuss 
in Section [V] the quality of the proposed wavefunction by 
comparing it with other variational wavefunctions. 



II. MODEL AND METHOD 

The Hamiltonian of the Hubbard-Holstein model at 
half-filling reads: 

H = -t ]T c\ a c ja + jJ2( Ul ~ + w oXX fli 

<i,j>,a i i 

+au>o y^,( n i - l)(°f + a 0i (!) 



where Ci a (c| CT ) and a ia {a\ a ) are respectively annihila- 
tion (creation) operators for tight-binding electrons with 
spin a and for optical phonons of frequency ujq on site 
i, t is the nearest-neighbor hopping amplitude, U the lo- 
cal Hubbard repulsion and a parameterizes the coupling 
between local displacements and electronic density fluc- 
tuations Uj — 1. Both interaction terms have been written 
in a particle-hole symmetric form which enforces the half- 
filling condition. We introduce an adiabaticity parame- 
ter defined as 7 — loq/D, measuring the ratio between 
the phonon and electron characteristic energies, and a 
dimensionless parameter A = 2a 2 u>a/D = 2a 2 ~f which 
measures the strength of the e-ph coupling. Here D is 
the half bandwidth, which depends on the hopping pa- 
rameter t once a particular lattice is chosen. We perform 
our calculations in an infinite-coordination Bethe lat- 
tice with semicircular density of states of half-bandwidth 
D = 2tM For this lattice, as in any infinite-coordination 
lattice, the expectation value of the Hamiltonian on the 
Gutzwiller state can be performed without further ap- 
proximations, because all the configurations with the 
same average occupation number equally contribute to 
the quantum averages. The same approach can be used 
also for finite-dimensional lattices, where it represents a 
further approximation, usually called Gutzwiller approx- 
imation. Another advantage of considering the infinite- 
coordination Bethe lattice is the possibility to compare 
with the exact results for this lattice, obtained through 
DMFTJ^±id£ 

We introduce the GPW as^L 

I*) = HiPi(xi)\* ), (2) 

where |^o) is a Slater determinant, that we take as the 
non-interacting Fermi sea in order to describe metallic 
states. Viipui) is a generalized Gutzwiller projection op- 
erator 

Vi{xi)= J-^ < M^)h)(^l- (3) 

i/=0,l,2 V ■ 

Here \vi){vi\ are the projectors associated to the different 
local electronic states, empty site (y — 0), singly occu- 
pied site (y = 1) and doubly occupied site (v = 2). Our 
wavefunction therefore associates to each local electronic 
state a normalized first quantization phonon wave func- 
tion (j) v {xi) depending on the displacement coordinate 
Xi = {fli + a\)/y/2, which has to be determined varia- 
tionally, exactly like the probabilities of each of the local 
states P„ (P^ ' are the same quantities for the uncorre- 
lated system, i.e., P (0) = (1 - n/2) 2 , P x (0) = n/(l - n/2), 
and — (n/2) 2 , n being the average density per site). 
By imposing the standard constraint s 1 ^ 34 

J dar t - <*o|^| a |*o) = 1> (4) 
J f&i^oklPiH^o) = n, (5) 
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and introducing the parameter d = (Pq + -P2V2 and 
the doping S = 1 — n, one gets P = d + | , Pi = 
1 — 2d and P 2 = d — | for the correlated occupation prob- 
abilities. Setting 5 = 0, the three amplitudes are con- 
trolled by a single variational parameter, d. The varia- 
tional energy per site is 



E = p »(M*))v + V2aLu [P Q (x)o 

v=0.1,2 

- 2| £ ||S| 2 + I (P 2 + P ) 
where h (x) = (o;o/2)(— d 2 4- x 2 )) and 

/» CO 

(0),= / d% <f>t(x)0<j> v (x) 



P2(xh] 



(6) 



(7) 



indicates the average over the phonon wave function 
4> u (x) of the operator O. The second line of Eq.© repre- 
sents the standard energy of the Gutzwiller approxima- 
tion for the Hubbard model with 



2|e| = t(* | E 4^l*o) 



(8) 



<i,j>,a 

the free electron kinetic energy and 



S = 



i/=0,l 



1 — J 2 



da; 4>* v+1 {x)<l> v {x) 



(9) 



the renormalization factor that accounts for the reduced 
electron mobility in the presence of e-ph and e-e interac- 
tions. Then we have to minimize with respect to d and 
the <j> v (x). The first minimization gives 



U + ((h (x)) + (h {x)) 2 - 2(h {x)) 1 ) 

.d\S\ 2 



+V2auj {{x) - (x) 2 )-2|e| 



i)d 



0, 



(10) 



while the second yields the following non-linear second- 
order differential equations 



4> a = h (x)4> + \/2auj x4>o - 



Pa 
ei_ 

Pi 

— 4>2 = ho{x)(j>2 



4>i = h (x)(f>i 



'S* 



V2auJoX(f>2 




where the e„'s are Lagrange multipliers enforcing the nor- 
malization conditions on the (/>„'s. 

Before discussing in some detail the solution of 
Eqs. (TTTj) - p3[) . we briefly analyze their structure. In the 
absence of e-ph coupling they reduce to three equivalent 
Schrodinger equations for a free harmonic oscillator cen- 
tered at each site independent on the electron occupancy. 
Therefore the phonon wavefunction is the standard gaus- 
sian, contributing ujq/2 to the on-site energy. Switching 
on the e-ph interaction has two major effects. The first 



is to induce a shift proportional to ±\/2a to the phonon 
wavefunctions associated to the v = 0, 2 electron states, 
which is the expected effect in the atomic limit with the 
chosen form of the e-ph interaction; the second is a non- 
local coupling between wavefunctions related to different 
charge configurations, which appears to be proportional 
to the average kinetic energy of the electrons. This is not 
surprising, since one expects the electron dynamics to af- 
fect the phononic properties driving them away from the 
atomic limit scenario. Nonetheless, it is readily seen that 
this coupling term is inversely proportional to the adi- 
abaticity ratio, that means, as it is well-known, that in 
the antiadiabatic regime, i.e. when phonons move faster 
than electrons, the interplay of phonon and electron dy- 
namics is negligible: the lattice rearranges itself almost 
instantaneously with respect to slowly moving electrons 
and a charge-dependent shift of the phonons is sufficient 
to capture the ground-state properties of the system. On 
the other hand the coupling term becomes dominant in 
the opposite limit, i.e. for 7 <^ 1; in this case each pro- 
jected wavefunction depends heavily on the other wave- 
functions and on the electron probability distribution, 
thus suggesting non trivial dependence on correlation ef- 
fects. 

At half-filling we have Po = P2 = d and Pi = 1 — 
2d, that imply (f>i(— x) = <fii(x) and <p2(—x) = (j>o(x) for 
the lowest-energy state. Therefore we are left with two 
variational equations for the lattice ground state coupled 
with Eq. (flTj|) . The latter can be recast in the more 
transparent form 

4d = 1 - u, (14) 
which depends from the important quantity 

U - 2 y/2 au>o {x)% — 2 (ho(x))i + 2 (ho(x))2 

U= 8| £0 ||(<^i>| 2 1 ' 

measuring the effective degree of electronic correla- 
tion as the ratio between renormalized U e ft (nu- 
merator) and phonon-renormalized kinetic energy s e ff 
(denominator)^. 

In spite of these simplifications, an analytical solution 
of our equations is still a difficult task. In the next section 
we present the general solution of Eqs. (fl~Tj) - (fT3)) . Here 
we briefly discuss some specific choices for the phonon 
wavefunction that allow for an analytical solution and 
reproduce popular variational approaches. 

a. Gaussian ansatz for the phonon wavefunctions. 
As discussed above, in the atomic limit the phonon wave- 
functions reduce to displaced harmonic oscillators whose 
displacement xq depends on the local electron occupation 
according to xq(v) = \/2a(y — 1). We can assume: 



|0) 



(16) 



where |0) is the ground state of an undisplaced harmonic 



oscillator, p 



-i(a 



f )/V2 is the conjugate momen- 



tum and f v is a parameter to be variationally deter- 
mined. Exploiting the particle-hole symmetry we can 
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impose f% = — fo = / and f\ — and obtain the varia- 
tional ground-state energy per-site: 

^ - M{l-2d)e- a2 f 2 |e | + 

d[[f-2aW(2-/)], (17) 



So 



that corresponds, within a constant, to the result of 
KR slave-bosons supplemented by the variational Lang- 
Firsov transformation (VLF) 28 i 29 . This is not surprising, 
since the Gutzwiller approach is known to be equivalent 
to the slave-boson mean-field approach and the VLF de- 
scribes the phonon wavefunctions as displaced Gaussians. 

The mean-field solution is then determined by min- 
imizing (fTT)) . that amounts to solve Eq. (fT4")) with u — 
exp(a 2 f 2 ) [U - 2a 2 uj / (2 - /)]/8|e | together with 2 ^ 



(18) 



b. Squeezed phonon state ansatz for the phonon 
wavefunctions. As widely discussed for the Holstein 
model— i 24 ' 26 i 27 , the harmonic (gaussian) ansatz is not 
expected to be accurate except for the so-called light po- 
laron case, that is realized when anharmonic fluctuations 
of the lattice are not so important, i.e. when A, 7 <C 1 or 
when 7 > 1, \U e fff*. Such anharmonic fluctuations can 
be partially captured by assuming that phonon wavefunc- 
tions are "squeezed" two-phonon coherent states (Varia- 
tional SQueezed, VSQ). The simplest choice that is usu- 
ally made is: 



1^) 



e -af v (a i -a) e -F{aa-a'< a f ) 



10), 



(19) 



where the e-ph induced displacement of the phonon field 
is still associated to the local charge state whereas the 
variational parameter F > controlling the squeezing of 
the wavefunctions is charge-independent, i.e., as noticed 
first in Ref. 24, it describes a kind of mean-field effect on 
the lattice due to the electron motion. The ground-state 
energy computed over these phonon wavefunctions reads: 



* - 



t- 2 ) - 8(2(1 - 2d)e 



-a 2 / 2 



£0 



d[f/-2aW(2-/)], (20) 

where we have introduced t 2 = exp(— AF). This is what 
one would obtain by applying the slave-boson mean field 
to the effective polaron model derived by averaging over 
squeezed phonon state a 23 ' 26 . The mean-field equations 
are modified as follows: 

-1 



1 



gN 



(1 + u) e 



_ r 2 2 

-a f t 



■ a 



2 f2 



f 2 (l~u 2 )e 



2 ?2 2 

-a f T 



-1/4 



(21) 



(22) 



with u = exp(a 2 / 2 r 2 ) [U-2a 2 u f(2- /)]/8 |e |. While / 
still accounts for the effective e-ph-induced displacement 
of the ions, the variational parameter r can be viewed 
as a measure of anharmonic fluctuations of the phonon 
wavefunctions. 



III. HALF-FILLING PHASE DIAGRAM. 

As we briefly discussed in the introduction, the param- 
agnetic phase of the half-filled Hubbard-Holstein model 
is an interesting test-field for our trial wavefunction since 
the interplay of the two local interaction mechanisms 
is more effective and at the same time more transpar- 
ent because both the interaction terms are able to drive 
metal-insulator transitions of different nature. From a 
technical point of view, imposing the particle-hole sym- 
metry simplifies the calculation, eliminating of the three 
self-consistency equations for the phonon wavefunctions. 
Finally, the half-filled phase diagram has been exten- 
sively studied by means of DMFT, providing us with 
an exact benchmark for the infinite coordination Bethe 
lattice^i 3 -^ 

In this section we mainly focus on the metal-insulator 
transitions driven by the two interaction terms. In the 
absence of e-ph coupling a transition from a metal to 
a Mott insulator (MI) occurs at a critical U c , whereas 
the half-filled Holstein model displays with increasing A 
a transition to a pair (bipolaron) insulator (BPI), the 
occurrence of such instabilities depending on the adia- 
baticity regime 3 ^. 

In the present framework, a key parameter to charac- 
terize the properties of the system is the effective corre- 
lation parameter u. From Eq. (| 14[) . u = 1 corresponds 
to d = 0, and u — —1 to d = 1/2. These are suffi- 
cient conditions for, respectively, the Mott and bipola- 
ronic transitions. In fact, in the present mean-field de- 
scription, the vanishing of the average double occupa- 
tion signals the transition to a MI, with one electron per 
site, whereas d — 1/2 corresponds to a system of elec- 
tron pairs stuck to lattice sites, i.e., a BPI. Once the 
existence of the insulating solution is proved, their ther- 
modynamic stability must however be checked by com- 
paring their ground-state energies (i.e., E (MI) = oj$/2 
and E a (BPI) = [uo + U- 2a 2 uj }/2) with that of other 
possible (metallic) solutions. To gain a first insight about 
the combined role of the two interaction terms, we notice 
that the effect of phonons on the correlation parameter 
u is twofold. On one hand the e-ph interaction reduces 
the numerator of u, but it also reduces the denominator 
normalizing the kinetic energy. This means that under 
different circumstances the two interacting mechanisms 
can either cooperate or compete in localizing the parti- 
cles. 

A full solution of the Gutzwiller equations can be eas- 
ily achieved numerically by expanding the phonon wave 
functions in the complete basis of the eigenstates \n) of 
ho, the harmonic oscillator centered around x = 0, with 
eigenvalues £ n = wo(n + 1/2) 



(23) 



n=0 



This expansion is completely general and it does not 
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introduce any further constraint or approximation. Ex- 
ploiting the particle-hole symmetry ci ' = (— l) n ci 2 ' ) anc 
plugging K23J) in the GPW equations ^TT])-((T3j) one gets: 



imposing u = 1 one gets for small A the explicit relations: 



(2 



(1) Jl)" 



-£ n )<W+4|e |(l-2d) cW C 
-V2. 



c (1) 



J 2 ) 



= 



which can be solved iteratively by keeping an arbitrary 
number of harmonic oscillator levels (in practice with 50- 
70 levels we get already very accurate results for the 
ground state); of course their solution depends on the 
value of the parameter d, which has to be determined 
self-consistently through Eqs. (fl4 |) .(ll5 |) at each iteration. 

We show in Fig. [T] the obtained phase diagram in the 
U — A plane, displaying three different phases, the two 
insulators MI and BPI and a correlated metal, which, as 
we will discuss in the next section, may show bipolaronic 
features before the metal-insulator transition. Since one 
expects the properties of the system to be strongly de- 
pendent on 7, as it happens in the absence of correla- 
tion, the phase diagram has been determined for three 
different values of the adiabaticity parameter, namely 
7 = 0.2,0.6,4. 

In the antiadiabatic limit, where the phonon wavefunc- 
tions are simply displaced harmonic oscillators, one finds 
that | (<?(>i |</> 2 ) | = exp(-a 2 /2) w 1 and u ~ [U-XD]/U C = 
U e f f/U c - The Holstein attraction becomes instantaneous 
and the system behaves essentially as a repulsive Hub- 
bard model as long as U e ff > and as an attractive 
one for negative U e ff. On the repulsive side we will 
have a Mott transition, while on the attractive side the 
normal phase displays a metal-insulator transition be- 
tween a metal and an insulating phase formed by in- 
coherent pairs, which is the antiadiabatic limit of the 
BPI 36 . Therefore transitions to insulating phases are sec- 
ond order and occurs symmetrically with respect to the 
U e ff = line, i.e. U M i = U C + XD and XbpiD = U C + U 
(thin dotted lines in Fig. [I}. 

Considering finite values of 7 we find that the main 
effect onto the Mott transition with respect to 7 — > 00 
is to change the slope of the transition line at small A, 
which remains second-order and is shifted to lower values 
of U with decreasing 7. Increasing A the metal-insulator 
line bends towards the bisector U = XD when 7 is de- 
creased. We notice that for small A the Mott transition 
line basically coincides with the prediction of the simple 
VLF approach. This is easily understood observing that, 
as long as the transition is of second order, it is defined 
by the condition u = 1. According to Eq. (|22|) . t 2 — > 1 
as u — > 1, which means that squeezing effects do not 
improve the VLF variational ansatz. Plugging the VLF 
expression for the phonon wavefunctions in Eq. (|15p and 



ZUJq 



XmiD = 



U 

2uj q 



7 « 1, (24) 



(U - U c ) 7 » 1, (25) 



in excellent agreement with previous DMFT phase 
diagram o 12 ' 13 ! 14 . 

Turning to the region of the parameter space domi- 
nated by the e-ph coupling, a more complex, strongly 
7— dependent picture emerges. As long as 7 » 1, the 
metal-BPI transition is found to be second-order, and 
the condition u = — 1 is always satisfied at the transi- 
tion. Since in this limit the (\> v 's are essentially displaced 
harmonic oscillators, one finds: 



XbpiD = U c 



e~- + U/U c 



(26) 



As expected, the presence of U obstacles the stabiliza- 
tion of the BPI and pushes the pair transition to higher 
values of A. The ability of U to compensate the phonon- 
mediated attraction is maximum for large 7, so that the 
same transition line moves to smaller A (yet larger than 
in the absence of U) as 7 is reduced. Decreasing 7 the 
scenario becomes more involved, and the order of the 
transition turns out to depend on both 7 and U . For 
0-5^7^1 there is a window of intermediate values of 
U in which the transition becomes of first order. There- 
fore the transition to the BPI is second-order at small 
values of U, turns to first order at a given U depending 
on the considered 7 and then, at larger U, the transition 
becomes second-order again (e.g. at 7 = 0.6 we obtain a 
first-order transition for U/U c between 0.5 and 1.1). The 
metallic region comprised between the MI and the BPI 
appears to close asymptotically on the line U e f / — with 
increasing interactions strength. Eventually, for 7 < 0.5 
the transition to the BPI is first order for all values of U. 

Some comments are in order about the evolution 
of the order of the transition to the BPI. For zero 
electron-electron interaction, DMFT shows a continuous 
transition^. Thus the first-order nature of the metal- 
BPI transition at zero and small U is due to the limita- 
tions of our variational approach. Indeed more restrictive 
variational choices, like the VLF and the VSQ predict a 
first-order transition for 7 < 71 (for the VSQ the limiting 
71 is only slightly smaller than for VLF). The full solu- 
tion of our GPW equations lowers the value of 71 , but a 
more sophisticated wavefunction, able to capture nonlo- 
cal character of the ground state of e-ph coupled system, 
is required in adiabatic regimes, when retardation effects 
between the motion of electrons and lattice distortion of 
the lattice become particularly relevant. 

On the other hand a change in the order of the tran- 
sition, from second to first one with increasing U, has 
been reported within exact DMFT at U/Uc — 0.5 for 
7 = 0.1—, and it was associated in Ref. [l5| to an abrupt 
change in the electronic configuration from the Mott state 
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where all the sites are singly occupied to the BPI where 
all the sites are either empty or doubly occupied. There- 
fore our GPW is able to correctly reproduce this highly 
non trivial evolution of the order of the transition as a 
function of U. We underline that more limited variational 
choices such as VLF and VSQ always find a first-order 
transition regardless the value of U. This result witnesses 
that the improvement brought by the GPW is particu- 
larly important in the presence of correlation effects. 

Our analysis shows a remarkable agreement between 
our phase diagram and previous DMFT results for 7 = 
0.2. Such an accuracy can only increase for larger 7, 
where we observe, e.g., a metallic phase which intrudes 
between the two insulators for large U and A. Unfor- 
tunately to our knowledge there are no detailed DMFT 
studies of the phase diagram for larger values of 7 to 
compare with our variational predictions. Recently, the 
phase diagram of the one-dimensional Hubbard-Holstein 
model at different values of the adiabaticity parameter 
has been obtained, and an intermediate metallic phase 
has been observed getting larger between the two insu- 
lating phases as 7 increases^, in qualitative agreement 
with our results. 



IV. QUASIPARTICLE WEIGHT, PHONON 
DISPLACEMENT AND BIPOLARON 
FORMATION 

Beside the metal-insulator transition associated to 
bipolaronic binding, a strong electron-phonon coupling 
can drive the formation of polaronic and bipolaronic 
states. A polaron is an electron strongly entangled 
with phonons, thus moving with a significantly enhanced 
effective mass. The residual interaction between po- 
larons is attractive, leading to the formation of pairs, 
called bipolarons. Strong Coulomb repulsion or large 
lattice fluctuations can avoid the metal-insulator tran- 
sition due to bipolaron localization, and they can lead 
to a poorly metallic state of bipolarons. It has been re- 
peatedly shown that the formation of polarons and bipo- 
larons in the absence of electron-electron interaction oc- 
curs through a continuous crossover. In the following 
we introduce and discuss two quantities that allow us to 
characterize the (bi)polaronic character of the correlated 
metallic phase from the point of view of both electronic 
and phononic observables. 

The metallic or insulating character of the variational 
solutions can be addressed by computing the quasiparti- 
cle renormalization factor, which corresponds to the in- 
verse effective mass m* in the Gutzwiller approach: 

Z = = 8d(l - 2d)|(0 1 |0 2 )| 2 = (1 - u 2 )K&|fc)| 2 , 

(27) 

that is zero when the system is a MI or a BPI and a 
nonvanishing quantity < 1 for metallic solutions. The 
effect of the e-ph coupling on the mass renormalization 
is not only comprised in the overlap (4>i\4>2), a s in the 



infinite-[/ limit 31 , but it is also hidden in the effective 
correlation parameter u. 

In the absence of U it is clear from Eq. (|27|) that the 
main correction to Z for weak e-ph coupling effects comes 
from the overlap (4>i\(f>2), which actually results in an en- 
hancement of the effective mass that is linear in A, as 
shown in the left panels of Fig. [2J however the rapid sup- 
pression of Z as the BPI is approached by increasing A 
mainly comes from the phonon-induced attractive inter- 
action, at least for 7 > 0.5 (top and middle panels of the 
left column in Fig. [5]). Approaching the adiabatic regime 
the relevance of the effective phonon-mediated attraction 
is reduced, with Z almost completely determined by the 
behavior of the overlap (<fii Ifa) as a function of A (bottom 
panel of the left column in Fig. [2]). Yet, the sudden de- 
crease of the overlap drives a sharp first-order transition 
to the BPI for 7 < 0.5. 

Moving to the correlated metallic phase with U = 
0.8U C , we find that, for all the considered values of 7, 
Z has a non monotonic behavior. For small coupling Z 
increases with A, it reaches a maximum before decreas- 
ing to eventually reach the BPI. Also in this case Eq.([27p 
allows to explain in a very intuitive way the observed 
increase of Z when the e-ph coupling is turned on in a 
correlated regim o 12 i 14 . In fact, as long as the effective 
interaction stays repulsive, the main effect of A is to par- 
tially screen the bare Hubbard repulsion, the effectiveness 
of such screening being strongly dependent on the adia- 
baticity regime 14 (cfr. the right panels of Fig. [2|); the 
electron motion, even if slowed down by e-ph coupling, 
can take advantage of the smaller energy cost associated 
to double occupation, and this results in a decrease of m* 
with respect to the A = value. As the phonon-induced 
attraction overcomes the Hubbard repulsion, the two ef- 
fect of phonons (localizing tendency to self-trap and to 
form local pairs) are cooperative, leading to a sharp first- 
order transition to the BPI for strong enough U and small 

71a. 

The above discussion clearly shows that Z is not the 
right quantity to identify the bipolaron crossover in the 
presence of competing interactions. While in the ab- 
sence of U bipolaron formation is usually associated to 
an abrupt enhancement of the effective mass, as shown 
in the left panels of Fig. [21 a more direct measure of po- 
laronic features in the groundstate has to be used when 
the effective mass results from the competition between 
attractive and repulsive interactions. More precisely, we 
need a quantity that measures the entanglement between 
the electron motion and the lattice distortion. A good 
candidate is the groundstate phonon distribution func- 
tion 

P(x) = (Vok)(2#o) = ^P„|<Ma;)| 2 , (28) 

\x) being an eigenstate of the displacement operators, 
and |V>o) the groundstate vector. P{x) measures the ef- 
fective (ground-state) displacement of the lattice. In fact, 
at least in the U — limit the development of polaronic 
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features is associated to the evolution from a monomodal 
shape of P(x) centered around an equilibrium displace- 
ment (x = at half-filling) to a bimodal distribution with 
maxima displaces of ±20, signaling the development of 
finite lattice distortions—. Focusing on the half-filling 
case, x = changes from a maximum in the monomodal 
distribution to a minimum in the bimodal one. Thus we 
can associate the appearance of bipolarons in our model 
to a change of sign in the second derivative of P{x) for 
x = 0, i.e., 



4rf 



dx 2 



x=0 



dx 



x=0 



2(1 -2d) (f> 1 (x = 0)- 



dx 2 



> 0. 



(29) 



x=0 



We notice that the evaluation of the condition (|2"9")l is 
particularly sensitive to the choice of the variational pro- 
jected wavefunctions. The evolution of P(x) for our 
Gutzwiller function is shown in Fig. [3] in comparison 
with the VSQ ansatz in a moderately adiabatic situa- 
tion 7 = 0.6. We notice that some care has to be taken 
in using the condition (|29p to pinpoint the bipolaron 
crossover in the correlated regime (right panels), since 
the effect of U on the P{x) can result in shoulders or 
even in a three-peak structure rather than the simple bi- 
modal shape. However our criterion captures the correct 
order of magnitude of the critical A for the onset of a 
bipolaronic groundstate. 

In the antiadiabatic limit 7 — > 00 no bipolaron 
crossover occurs before the system turns insulating by 
forming incoherent local pairs: Replacing <p v in Eq. (|29|) 
with the proper displaced harmonic oscillators, one finds, 
in agreement with the DMFT of Ref. HH, that bipo- 
laron formation occurs at the transition line only when 
a 2 > 1/4, which implies a very large A when 7 is large. 
However, at finite values of the adiabaticity parameter, 
the condition (|29p can be realized, since the lattice dis- 
tortions induced by the e-ph coupling, which are roughly 
proportional to 1/^/7, can push the system to polaron 
formation before the BPI instability takes place. Actu- 
ally, for all the considered values of 7 we found bipo- 
laron formation before the metal-BPI transition line, in 
a region whose size depends again on both correlation 
and adiabaticity parameter. On the weakly correlated 
side, we notice (see Fig. [T]) that, starting from small 7, 
the distance between the bipolaronic crossover and the 
metal-insulator transition first increases before decreas- 
ing in the antiadiabatic regime. This is in qualitative 
agreement with DMFT^, even if for 7 = 4 the GPW still 
predicts the bipolaronic crossover before the BPI transi- 
tion, in contrast with DMFT, where for the same value 
of 7 the opposite occurs. As shown in Fig. [TJ at 7 = 4, 
the presence of U slightly enlarge the region where the 
ground state displays polaronic features. 



V. DISCUSSION OF THE RESULTS. 

In this section we discuss the behavior of the phonon 
wavefunctions, that represent the main novelty intro- 
duced by our Gutzwiller method with respect to previous 
schemes where the functional form was assumed a priori. 

A comparison with VLF and VSQ schemes helps us to 
benchmark our method. Due to the enlarged variational 
space, the ground-state energy of the GPW is always 
smaller than that given in Eqs. (fl7|) and ([20]) . The gaus- 
sian ansatz for the ^„'s (VLF) results in a significantly 
higher energy, especially in the intermediate to strong 
coupling regime, whereas the VSQ energy stays very close 
to the general GPW one, even though the GPW becomes 
visibly more accurate than VSQ in the presence of size- 
able correlations, as shown in Fig. 2] However, even 
when GPW provides only a slight improvement on the 
variational energy with respect to VSQ, it gives a more 
accurate description of the physical properties of the lat- 
tice; in fact, the projected phonon wavefunctions of GPW 
can provide a better estimates of the different energy con- 
tributions, i.e. the potential energy gain and the kinetic 
energy loss due to lattice distortion, whose balance re- 
sults in the overall ground-state energy. 

To get a flavor of the improvement introduced by the 
wavefunction Eq. j2]) we plot in Fig. [5] 4>2 and tj>i as 
obtained by means of our general equations, compar- 
ing them with the simpler variational ansatzs, VLF and 
VSQ. In the absence of Coulomb repulsion, the squeezed 
states stay relatively close to the self-consistent solution, 
whereas VLF underestimates the effective lattice distor- 
tion and therefore the potential energy gain. On the 
other hand, in the strongly-correlated regime the GPW 
functions are rather different from the other methods and 
underline the non-trivial way in which anharmonic fluc- 
tuations of the lattice states are influenced by the corre- 
lated electron motion. Even if significantly more accurate 
than VLF, the VSQ underestimates the anomalous lat- 
tice fluctuations that we find in GPW; since the overlap 
of the phonon wavefunctions related to different charge 
configurations is essential in capturing the way e-ph cou- 
pling rcnormalizes the effective hopping, this means that 
VSQ overestimates the kinetic energy loss. 

A deeper understanding of the differences between the 
standard variational ansatzs and the self-consistent so- 
lution for the phonon wavefunctions can be traced out 
by looking at the degree of effective displacement, given 
by {x) v , and of the related quantum fluctuations, that 
can be estimated by Ax 2 — ([x v — (x v )] 2 ). These two 
quantities are, roughly speaking, related respectively to 
the potential energy gain and to the kinetic energy loss 
induced by the e-ph coupling, and their inspection will 
help us to identify the way in which the GPW improves 
with respect to the other approaches. It is readily found 
that (x) v = ±\/2af for v = 0, 2 and {x)\ = for 
both VLF and VSQ ansatz, with / given respectively 
by Eqs. (JTSJ) and (|2T]l. whereas Axl = 1/2 in the har- 
monic approximation and Ax 2 = l/2r 2 for any v in 
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the squeezed phonon state approximation. In Fig. [6] 
we compare the effective displacement and Ax 2 obtained 
within the different approaches as a function of A in the 
presence of a sizeable U, namely U — 0.8U C . By look- 
ing at (x)„, we find that both harmonic and squeezed 
ansatzs reproduce the weakly interacting regime quite 
accurately, while VLF deviates in the intermediate to 
strongly-coupled regime, falling suddenly onto the atomic 
solution (which in the present mean-field approach actu- 
ally translates in the BPI). The origin of such a discrep- 
ancy is the relevance of phonon quantum fluctuations 
that are crucial in the crossover region and are poorly 
described within VLF. However within GPW such effect 
is strongly dependent on the charge state associated to 
the different wavefunctions: the phonon wavefunctions 
projected onto empty and doubly-occupied sites display 
a maximum in Ax^n) an d then tend to the atomic so- 
lutions as the BPI is approached, whereas Ax\ increases 
monotonically. Therefore, as anticipated, the simple two- 
phonon coherent state Eq. (|19p reproduces only a kind 
of average squeezing effect; to stress this point we show 
an averaged Ax 2 , namely Ax 2 = ^2 PvAx 2 ,, that indeed 
stays very close to the VSQ solution. 

In conclusion, the improved accuracy of our variational 
phonon wavefunctions proves to be relevant especially for 
intermediate-to-strong e-ph coupling and in the presence 
of sizeable correlation, in particular close to the BPI in- 
stability, when quantum fluctuations of the phonons play 
an important role. Indeed, by looking at Fig. [6l one 
immediately realizes that the evaluation of the phonon 
distribution function strongly relies on the quality of the 
projected wavefunctions, and the fulfillment of the con- 
dition (|29f can lead, as it actually does, to completely 
different bipolaron crossover lines for different variational 
ansatzs when considering adiabatic regimes. To be more 
explicit, we find that the simple variational Lang-Firsov 
(or harmonic) ansatz never predicts the onset of a bipo- 
laronic groundstate for any U at 7 < 1, being P(x) al- 
ways monomodal (results not shown) . On the other hand 
VSQ wavefunctions reproduce a bimodal P(x), displayed 
in the bottom panels of Fig. [3J at least at U = 0, but at 
larger A with respect to GPW; however they completely 
miss the complex structures that P{x) develops in the 
presence of a sizeable U, and they do not reproduce any 
bipolaron crossover in the correlated regime. 



VI. CONCLUSIONS 

In this work we have extended the analysis of a new 
type of Gutzwiller variational wavefunction introduced 
to describe correlated electron systems in the presence 
of e-ph coupling. In a previous paper— we benchmarked 
the wavefunction in the infinite-?/ limit by comparing it 
to exact DMFT calculations, showing that polaron for- 
mation occurs smoothly for any value of the adiabaticity 
parameter 7 and electron density n, as opposed to the 
adiabatic first-order transition found in standard varia- 



tional approache d i 29 . In the infinite-C/ limit the effect 
of correlation on polaronic physics can be analyzed by 
tuning electron filling, and we found that polaron forma- 
tion is inhibited and pushed to higher, though finite, A 
when moving from low-density to half-filled system. 

Relaxing the assumption of infinite repulsion between 
electrons, in this paper we turned our attention to the 
half-filling regime and considered, together with bipo- 
laron formation, the competition between U and the 
phonon-mediated bipolaron instability, ruled out by an 
infinite U. We found, in excellent agreement with pre- 
vious work a 14 ' 28 ! 29 , that the Mott metal-insulator transi- 
tion is always robust with respect to polaron formation 
and that screening of the bare repulsion due to the cou- 
pling with the lattice is less and less effective as adiabatic 
regimes are attained. On the other hand we observed 
that both bipolaron crossover and metal-BPI transition 
depend not only on the adiabaticity parameter but also 
on the strength of Hubbard repulsion. Two localizing 
mechanisms are involved at intermcdiate-to-strong e-ph 
coupling, one resulting in self-trapping of electrons cou- 
pled with Holstein phonons, and the other related to the 
phonon-mediated e-e attraction. If on one hand correla- 
tion pushes to larger values of A both the onset of bipola- 
ronic features and the metal-BPI transition, it influences 
the two processes in a different way, which depend on 
the adiabaticity regime, giving rise to different scenar- 
ios. For example, when 7 > 1, U is more effective in 
screening the attractive interaction than in preventing 
the bipolaron crossover—. On the other hand at small 7 
the bipolaron crossover is inhibited by U faster than bipo- 
laronic transition, and polaronic features are observed in 
a narrow region that rapidly shrinks with increasing U; 
at the same time a change in the order of the transition 
occurs depending on the strength of correlation. 

A critical comparison with standard variational ap- 
proaches to Holstein phonons shows that the proposed 
wavefunction can be viewed as a proper generalization 
apt to include correlation effects, underlying the impor- 
tance of a reliable variational description of anharmonic 
fluctuations in the presence of U. Our wavefunction im- 
proves on previous variational approaches already in the 
absence of electron correlations (VLF and VSQ can in- 
deed be obtained as restricted solutions of our method 
limiting the variational freedom), but it still has some 
limitations in the deep adiabatic limit due to the nature 
of our phonon wavefunctions, that are associated to the 
different local electronic states. For the same reasons, 
our method becomes considerably more accurate in the 
presence of strong correlations, that make the electronic 
physics more local. 
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FIG. 1: Phase diagram of the half-filled Hubbard-Holstein 
model in the U — X plane as obtained in GPW for three dif- 
ferent values of 7. Solid lines represent the metal-insulators 
transitions, dashed lines are the bipolaron formation lines ob- 
tained by imposing condition Eq. (|29p as discussed in Section 
IIVI (the inset shows the narrow region where a bipolaronic 
metal is found for 7 = 0.2). Thin dotted lines are the antia- 
diabatic (7 — > 00) transition lines symmetric to the U = \D 
line, also displayed. The vertical arrows in the middle panel 
mark the boundary of the region in which the transition to 
the BPI is of first order. 
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FIG. 2: (Color online) Typical evolution of Z as a function 
of A with and without Hubbard repulsion. Parameters are 
7 = 4, 0.6, 0.2 from top to bottom and U = 0, U = 0.8U C 
in the left and right panel respectively. Open circles are the 
phonon contribution to the hopping renormalization coming 
from |<0i|0 2 }! 2 . 
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FIG. 3: GPW results (top panels) for P(x) as a function of 
e-ph coupling constant at 7 = 0.6 in the presence of Hubbard 
repulsion (right panel, U/U c = 0.8 and A ranging from 2.8 to 
3.12) and with U = (left panel, A ranging from 1 to 1.3). 
The bottom panels show the lattice distribution function as 
obtained with VSQ for the same parameters. 




-0.52 1 1 1 1 1 - 

0.2 0.4 0.6 0.8 1 
X 



-0.015 



-0.02 




- VLF 




-0.025 




VSQ 




-0.03 


— 


GPW 




-0.035 


■ " s *f-^-r~" — — 






-0.04 








-0.045 




''x 




-0.05 


U/U c =0.8 


\ 




-0.055 








-0.06 









1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 
X 



FIG. 4: (Color online) Comparison of the variational on-site 
energies obtained for the paramagnetic solution in the absence 
of Hubbard repulsion (upper panel) and close to the critical 
value U c = 8|eo| for the purely electronic Mott transition 
(lower panel) for 7 = 0.2. 



13 




FIG. 5: (Color online) Comparison of the projected phonon 
wavefunctions for 7 = 0.2 as obtained with GPW and with 
the variational ansatzs discussed in the text (top panel is U = 
0,A = 1; bottom panel is U = 0.8U C ,X = 2.8). The thin 
dotted line is the corresponding atomic solution. 




FIG. 6: (Color online) Comparison of the effective displace- 
ment and width of projected phonon wavefunctions in the 
metallic solution before the BPI is reached as obtained with 
GPW and with the variational ansatzs discussed in the text 
(parameters are 7 = 0.6 and U = 0.8U C ). The thin dotted 
line is the atomic solution. 



